#include <fstream>
#include <iostream>
#include <iomanip>
#include "TFile.h"
#include "TTree.h"
#include "TLine.h"
#include "TSystem.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TCanvas.h"
#include "TEllipse.h"
#include "TClonesArray.h"

using namespace std;

void Pos(){

  int ret = gSystem->Load( "../Objects/libObjectClass.so" );
  if( ret!=0 ){
    cout << "Failed to Load Library: ret = " << ret << endl;
    return;
  }

  const int NEV=512;

  float Er = 1;
  int       Ng = 40; // N generated photon/keV
  float EffTPB = 0.1;
  float  totQE = 0.25;

  float eff=Er*Ng*EffTPB*totQE;

  TFile *fread=new TFile("/scratch/shinichi/OutputFiles/100kg_scan_WithRef_UpperPMT.root");

  TTree *tree=(TTree*)fread->Get("tree");

  const int nPMT=38;

  TClonesArray *JPMTa[nPMT];
  int nJPMThits[nPMT];
  int nPhoton;
  for(int iPMT=0;iPMT<nPMT;iPMT++){
    JPMTa[iPMT]=new TClonesArray("PMT");
    tree->GetBranch(Form("JPMT%db",iPMT))->SetAutoDelete(kFALSE);
    tree->SetBranchAddress(Form("JPMT%db",iPMT),&JPMTa[iPMT]);
    tree->SetBranchAddress(Form("nJPMT%dhits",iPMT),&nJPMThits[iPMT]);
  }
  tree->SetBranchAddress("nPhoton",&nPhoton);

  TH1F *hJPMTx=new TH1F("hJPMTx","hJPMTx",100,-150,150);
  TH1F *hJPMTy=new TH1F("hJPMTy","hJPMTy",100,-150,150);
  TH1F *hJPMTz=new TH1F("hJPMTz","hJPMTz",100,-650,-550);

  TH2F *hPosEff_All=new TH2F("hPosEff_All","hPosEff_All",25,0,250,75,-900,-150);
  TH3F *hPosEff3D=new TH3F("hPosEff3D","hPosEff3D",25,0,250,25,0,250,75,-900,-150);
  TH2F *hPosEff_J[nPMT];
  //TH3F *hPosEff3D[nPMT];
  TH1F *hMeanJ[nPMT];
  for(int iPMT=0;iPMT<nPMT;iPMT++){
    hPosEff_J[iPMT]=new TH2F(Form("hPosEff_J%d",iPMT),Form("hPosEff_J%d",iPMT),
			     25,0,250,75,-900,-150);
    hMeanJ[iPMT]=new TH1F(Form("hMeanJ%d",iPMT),Form("hMeanJ%d",iPMT),
			  1000,0,1*Er);
  }
  TH1F *hMeanAll=new TH1F("hMeanAll","hMeanAll",1000,0,1*Er);

  int nev = tree->GetEntries();
  cout<<"NEV = "<<nev<<endl;

  for(int iev=0;iev<nev;iev++){

    if(iev%1000==0) cout<<"--- "<<std::setw(5)<<iev<<" Event ---"<<endl;

    for(int iPMT=0;iPMT<nPMT;iPMT++){
      JPMTa[iPMT]->Clear();
    }

    tree->GetEntry(iev);

    for(int iPMT=0;iPMT<nPMT;iPMT++){

      int nhit = nJPMThits[iPMT];
      for(int ihit=0;ihit<nhit;ihit++){
	PMT *jpmt = (PMT*)JPMTa[iPMT]->At(ihit);
	hJPMTx->Fill( jpmt->xpos );
	hJPMTy->Fill( jpmt->ypos );
	hJPMTz->Fill( jpmt->zpos );
      }

    }

    int ix = iev % 21;
    int iy = iev / 21 % 21;
    int iz = iev / 21 / 21 % 61;
    float xx = 0.+ix*10.;
    float yy = 0.+iy*10.;
    float zz = -805.+iz*10.;
    float R = sqrt( xx*xx+yy*yy );
    //----- Fiducial Cut
    if( R>200. ) continue;

    int SumPMT=0;
    for(int iPMT=0;iPMT<nPMT;iPMT++){
      hMeanJ[iPMT]   ->Fill( nJPMThits[iPMT]/(float)nPhoton );
      SumPMT += nJPMThits[iPMT];
    }
    hMeanAll ->Fill( SumPMT/(float)nPhoton );
    hPosEff3D->Fill( xx,yy,zz,SumPMT/(float)nPhoton );

    //----- Cross section
    if( yy!=0. ) continue;

    SumPMT=0;
    for(int iPMT=0;iPMT<nPMT;iPMT++){
      hPosEff_J[iPMT]->Fill( xx,zz,nJPMThits[iPMT]/(float)nPhoton );
      SumPMT += nJPMThits[iPMT];
    }
    hPosEff_All->Fill( xx,zz,SumPMT/(float)nPhoton*20/38.);

  }


  gStyle->SetOptStat(000000);

  cout<<"  All Eff(Mean) = "<<100*hMeanAll->GetMean()<<" %"<<endl;
  for(int iPMT=0;iPMT<nPMT;iPMT++){
    cout<<Form("JPMT%d Eff(Mean) = ",iPMT)
	<<100*hMeanJ[iPMT]->GetMean()<<" %"<<endl;
  }

  TCanvas *c1=new TCanvas("c1","c1",800,600);
  TH1 *frame1=gPad->DrawFrame(0,0,0.03*Er,3200);
  frame1->SetXTitle("Acceptance");
  frame1->SetYTitle("NEV");
  hMeanAll->SetLineWidth(2.5);
  hMeanAll->Draw("same");

  /*
  for(int iPMT=0;iPMT<nPMT;iPMT++){
    hMeanJ[iPMT]->SetLineColor(iPMT+1);
    hMeanJ[iPMT]->SetLineWidth(2.5);
    hMeanJ[iPMT]->Draw("same");
  }
  */

  TCanvas *c2=new TCanvas("c2","c2",1200,1000);
  c2->Divide(5,4);  
  float max = 0.20;

  for(int iPMT=0;iPMT<nPMT;iPMT++){
    c2->cd(iPMT+1);
    hPosEff_J[iPMT]->SetMaximum( max );
    hPosEff_J[iPMT]->SetXTitle("r (mm)");
    hPosEff_J[iPMT]->SetYTitle("Z (mm)");
    hPosEff_J[iPMT]->Draw("colz");
  }

  TCanvas *c3=new TCanvas("c3","c3",800,600);
  hPosEff3D->SetLineColor(4);
  hPosEff3D->Draw("LEGO2");

  TCanvas *c4=new TCanvas("c4","c4",450,1500);
  hPosEff_All->SetMaximum( max );
  hPosEff_All->SetXTitle("r (mm)");
  hPosEff_All->SetYTitle("Z (mm)");
  hPosEff_All->Draw("colz");

}
